Quantum oscillations of the specific heat in d-wave superconductors with loop current 

order 
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We report numerical results of quantum oscillations of the specific heat in the vortex state of a 
d x 2 _ H 2-wave superconductor in the presence of loop current ordev^, which gives rise to Fermi pockets 
coexisting with nodal d x 2_ y 2-w&ve superconductivity. Within a lattice tight-binding model, we find 
that in an intermediate temperature range, the oscillations seem to approximately follow Onsager 
relation with an effective charge comparable to the electric charge. However, the quasiparticle 
spectrum does not resemble Landau levels. In order to understand the origin of the oscillations, 
we also perform Franz- Tesanovic transformation in the presence of the loop order and find that in 
addition to scalar and Berry potentials^, one component of the gauge invariant superfluid velocity 
couples to the low lying Dirac particles as a component of a vector potential. The magnetic field 
associated with this vector potential vanishes on average but is highly non-uniform in the magnetic 
unit cell. We also compare the results with the model without the loop order but with Zeeman-like 
coupling which also induces Fermi pockets in the superconducting state. 



I. INTRODUCTION 

Coexistence of d-wave superconductivity and Fermi 
pockets in underdoped high temperature cuprate super- 
conductors has been suggested by recent quantum oscil- 
lation experiments^ - — . Whether the Fermi pockets are 
electron-like or hole-like, and whether there is one or 
more than one pocket, is still under intense debate. In 
the present work, we focus on quantum oscillations of 
the specific heat, measurements of which have been pre- 
sented in Ref4. With the application of a magnetic field 
H, the non-oscillatory component of the Sommerfeld co- 
efficient 7(if) of the ultrapure YB2CU3O6.56 exhibits \[H 
behavior, which is consistent with the d-wave supercon- 
ductivity in the vortex state. Remarkably, this field de- 
pendence persists well into the resistive state. In ad- 
dition, there are several signatures of the existence of 
Fermi pockets. First, the zero field Sommerfeld coeffi- 
cient, 7(0) ~ 1.9mJ/molK 2 , is finite, indicating finite 
density of states at zero energy in zero field. Note that 
the low energy quasiparticles (QP's) of d K 2_ a 2-wave su- 
perconductors are characterized by the linear Dirac-like 
dispersion near four nodal points. This results in linearly 
vanishing density of states at zero energy. Although the 
finite density of states may in principle be induced by 
impurity disorder, the YBCO samples under study 8 are 
believed to be too pure to account for the measured value 
of 7(0). The high purity is consistent with the observa- 
tion of the quantum oscillations as well as the extracted 
values of the Dingle temperature. Therefore, the physi- 
cal origin of the nonzero 7(0) is most likely intrinsic to 
this system. Second, the oscillatory component of j(H) 
exhibits quantum oscillations in high magnetic fields, pe- 
riodic in 1/H, which can be well fitted^ - — by Lifshitz- 
Kosevich (LK) formula. 

Such phenomenology is quite remarkable. On the one 
hand, the quantum oscillations appear to be due to Lan- 



dau level quantization of the electron orbits, and indi- 
cate the existence of Fermi pockets, while on the other, 
7(i?) ~ \IH is a signature of d x 2_ y 2-wave superconduct- 
ing gap and the vortex state. 

The origin of Fermi pockets in the superconducting 
state has been under debate. One possible scenario is 
that the Fermi pockets arise from the one-dimensional 
CuO chains hybridized into the BaO layers 8 -. If not 
gapped by a proximity effect down to the lowest temper- 
atures (~1K) at which 7(0) was extracted, such Fermi 
pockets could result in a finite 7(0), as well as quantum 
oscillations of j(H). This would account for the main 
experimental features. 

In this paper we critically examine another scenario in 
which loop current order induces Fermi pockets in the 
d x 2_ y 2 superconductor— In such an ordered phase, 
charge currents circulate within each unit cell (as shown 
in FigQ}, breaking time reversal symmetry and inversion 
symmetry, but not their product or the discrete trans- 
lational symmetry of the lattice. Two Fermi pockets of 
Bogoliubov quasiparticles, one electron-like and one hole- 
like, are formed, giving rise to nonzero density of states 
at zero energy^. This may account for the finite 7(0) 
and the question is whether it can also cause quantum 
oscillations of "f{H) in high fields. In a superconductor, 
the Bogoliubov QP's are linear combinations of electrons 
and holes, and therefore do not carry definite charge. 
On the other hand, the QP's in Fermi liquids do carry 
definite charge^. Therefore, it is not a priori obvious 
whether there are any quantum oscillations at all, and if 
yes, whether the oscillations obey Onsager relation as in 
Fermi liquids. In this paper, we show that the effective 
magnetic field experienced by the Dirac quasiparticles in 
the loop order state vanishes on average and does not 
lead to Landau quantization. As such any oscillations 
do not follow the detailed LK phenomenology. Neverthe- 
less, the effective magnetic field experienced by the Dirac 
quasiparticles is highly non-uniform and in an interme- 



diate temperature range the quantum oscillations of the 
specific heat appear to approximately obey Onsager rela- 
tion, with an effective charge comparable to the electric 
charge. 

We investigate the oscillations of the specific heat in 
both the tight-binding lattice formulation and in the con- 
tinuum formulation. We assume that the vortices form 
a square Abrikosov lattice. In the tight-binding lattice 
formulation, we take the vortices to sit inside the pla- 
quettes of the two-dimensional Cu02-like plane. In each 
magnetic unit cell, there are two singly quantized vor- 
tices with flux hc/2e. Since the vortices are placed at 
the centers of the plaquettes, the vortex lattice has to be 
commensurate with the underlying tight-binding lattice. 
This prevents us from sweeping the magnetic field contin- 
uously. Instead, in this case, we sweep the (Bogoliubov) 
Fermi pocket area by varying the overall magnitude of 
the loop current order in fixed magnetic fields, and in- 
vestigate the dependence of the density of states and the 
specific heat on the Fermi pocket area. The results for 
the density of states are shown in Fig. (J9j) , where we also 
show that they clearly differ from the density of states of 
Landau quantized anisotropic Dirac fermions. Neverthe- 
less, as shown in Figs.(@][7]), we find that in an intermedi- 
ate, magnetic field dependent, temperature window, the 
specific heat exhibits oscillations as a function of Fermi 
pocket area for the four values of the magnetic field stud- 
ied, ranging from 7.7T to 35. 6T. At the same time, the 
non-oscillatory component of j(H) does not follow \fll 
behavior (see Fig.©). 

To further understand the origin of this effect, we com- 
plement the tight-binding calculations with an approxi- 
mate continuum formulation. To this end, we linearize 
the Hamiltonian in the vicinity of the four nodal points, 
perform the Franz- Tesanovic transformation^, calculate 
the quasiparticle spectrum numerically using plane- wave 
diagonalization and calculate the specific heat. While 
we are well aware of the subtleties with the large gauge 
invariancei^ we are merely interested in the overall qual- 
itative aspects of the results and their dependence on the 
strength of loop order and magnetic field. We find that 
the result obtained using this second method is consistent 
with the one obtained in the tight-binding lattice formu- 
lation. The second method offers an additional advantage 
in that the external magnetic field can be changed con- 
tinuously; the resulting oscillations of specific heat are 
shown in Figs. qi2fl!3|) . 

Finally, we compare these results with the results ob- 
tained by varying the Zeeman energy but without loop 
current order. The Zeeman term shifts all four nodal 
points, resulting in four Fermi pockets. In this case, the 
oscillations do not obey Onsager relation at all. 

Our paper is organized as follows. In Sec. II, we set up 
both the lattice and the continuum Hamiltonians, and 
calculate the zero field spectrum. In Sec. Ill, we present 
the numerical results for the quantum oscillations of spe- 
cific heat as a function of the loop current order and 
Zeeman energy, and the density of states in the lattice 



formulation. In Sec. IV, we present the numerical results 
for the oscillations as a function of the loop current or- 
der and the magnetic field in the linearized problem. In 
Sec.V, we discuss our results. 



II. FORMALISM: BDG HAMILTONIAN AND 
SINGULAR GAUGE TRANSFORMATION WITH 
LOOP CURRENT ORDER 

A. Lattice formulation 

We model the Cu02 plane in YBCO as a tight-binding 
lattice with lattice constant a, which may be set to I 
for convenience. (When converting to real units, we use 
a = 0.38nm.) When an external magnetic field H in 
the range H c \ < H < H C 2 is applied, the ci-wave super- 
conductor enters vortex state and the vortices form an 
Abrikosov lattice. We assume that a square vortex lat- 
tice is formed with magnetic unit cell £b x £b, where the 
magnetic length £b is define d through the flux quantum 
<fio = hc/e as £b = yj4>o/H. In each magnetic unit cell, 
there are two singly quantized vortices, each of which car- 
rying flux hc/2e. Our starting point is the Hamiltonian 
with nearest neighbor hopping, d- wave pairing and loop 
current order on the underlying tight-binding lattice in 
the presence of a magnetic field, 

H = H + Uj (I) 

where 14 

Ho = -t (e" M "'4 CT c r / (7 +/i.c.) 

{rr') a 

+ ( A rr< + 4, f C^) + h.c) (2) 

<rr') 

and the Hamiltonian for loop current order is 10 

Hj = Y(-iJrr'ci a c r , a + h.c). (3) 

rr'a 

In Eq.Q, the sums are over nearest neighbors (rr'), and 
a denotes the spin. The d-wave signature is encoded 
in i]s — +(— ) if <5||x(y). In the symmetric gauge, the 
magnetic flux $ through an elementary plaquette enters 
the Peierls factor via A rr+ x = — 7ry<f>/0o and A rr+ y = 
ttx^/c/jq. The d-wave pairing field in the vortex lattice 
is A rr / = r/ r ^ r r Aoe 1 ®"' , where the Ansatz for the pair 
phases isM. 

#(r) _|_ p #(r') 
~ | e #(r) +e #(r')|' W 

where V x V0(r) = 2ttz £\ 8(r - r*) and V ■ V0(r) = 
where Yi denotes the vortex positions. In Eq.©, the con- 
nectivity of the loop current network J rr / is determined 



2 



according to FigQ] and in zero field all the nonzero cur- 
rents have the same magnitude J 10 i 15 . In a finite mag- 
netic field, we have, explicitly, 



n, = -i 



"EG 



c rtT C r+XCT 



+xcr L r+y(T 



where A r+ x, r +y = —tt(x + y + 1)$ 
particle-hole transformation, 



if' 



■rt 

J 



Then the diagonalization of the H 
lent to the solution of the Bogoliul 
equation Hip r = Erp r where the lat 



n = 



A; 



Ar 

-E* + J* 



(8) 



with n being the chemical potential. The BdG Hamil- 
tonian acts on the two component Nambu spinor tj) T — 
[u r , v r ] T , and £ , J and A are defined through their action 
on a lattice function / r as 



<5=±x,±y 

J r fr = ~lJ Y, e_lArr+5 /r+^ 
5=x,— y,y— x 

(5=±x,±y 



(9) 
(10) 
(11) 



The Hamiltonian is invariant under discrete trans- 
lations followed by a gauge transformation (magnetic 
translations). As shown in RefJ^, it can be transformed 
into a periodic Hamiltonian by a singular gauge transfor- 
mation 



U 









-i<t>hM 



(12) 



where e (r) and </>/i(r) satisfy e (r) + </>h(r) = 0(r). The 
vortices are divided into two groups A and -B, and each 
magnetic unit cell contains one A and one B vortex, as 
shown in Fig[2] Then two phase fields ^(r) and </>s(r) 
are identified with <p e (r) and (fih( r )i respectively. If we 
choose 4>a{y) = 0s(r) = <fi(r)/2, then the transformation 
becomes IA = exp ^a^4>{r). Connecting pairs of vortices 
in one magnetic unit cell by a branch cut as shown in 
Fig 13 we have 



(13) 



where as discussed in detail in Ref J^, the Z 2 field 2:2, rr ' = 
1 on each bond except the ones crossing the branch cut 





FIG. 1: (Upper) The connectivity of the loop current net- 
work. (Lower) The dispersion of d-wave superconductor with 
loop current order in the first Brillouin zone. For clarity, the 
anisotropy is set to 1, and the strength of loop current order 
J=0.5. The dashed line indicates the zero energy. 
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FIG. 2: Magnetic unit cell Ib x £b containing A and B vortex 
joined by a branch-cut with Ib = 6a. 



where Z2. rr > = — 1. 
H = U^HU is 



Then the transformed Hamiltonian 



H = (T 3 (£ r - M) + ^lAr + JrX 



(14) 
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where the transformed lattice operators satisfy 

£ T ip r = -t J2 z ^-+ s x e <<73Vrr+ ^ r +5, (15) 
<5=±x,±y 

J r ^ r = -U ^2,rr + iXe lff3V -+ s i +j ,(16) 
<5=x,-y,y-x 

Arl/'r = A ^ Z 2,rr+<5 X rj S tp r +S, (17) 
<5=±i,±y 

cr's are Pauli matrices and 1 is the identity matrix, and 



|l + e «(^(r)-^(r'))|' 



(18) 



The resulting Hamiltonian is invariant under magnetic 
translations by £b in both directions, so it can be diago- 
nalized in the Bloch basis. The transformed Hamiltonian 
H{k) = e - lkr "He ikr becomes 

W(k) - CT 3 (£r(k) - M ) + ^lAr(k) + Jr(k)l (19) 

where 

£.(k)Vr = -t ^rr+ 5 X e^^+V^r+S, (20) 

<5=±x,±y 

Jr(k)V r = -U ]T Z ^+ S x e^^V^^l) 

5=*,— y,y-x 

A r (k)?/> r = A ^ z 2:rr+(5 x r/se^ipr+s. (22) 

<5=±x,±y 

B. Zero-field spectrum 

In the absence of a magnetic field, the phase factors 
e iA rr j an( j e iV rr / b ecome i f anc [ the Hamiltonian can be 

easily diagonalized, with eigenenergies 



-Ek = ±y + A 2 . + 2 J[sin k x — sin k y + sm(k y — k x Jft3) 

where = 2t(cosk x + cosk y ) — n and Ak = 2A (cos fc x — 
cos k y ). In the case with J = 0, the four nodes of the 
spectrum are located at (±fcu,±fe.D) where 



kn = arccos ( — ). 

Mr 



(24) 



In the vicinity of each node, the dispersion can be lin- 
earized 



Ei 



= y/v 2 F Sk 2 ± + V 2 A Skf 



(25) 



where 5k±{8kn) is the displacement of the momentum 
from a node in the direction perpendicular (parallel) to 
the Fermi surface, and the velocities are 



r /? = 2v^/l-(f) 2 t, 



r A = 2V2 x /l - (J) 2 A , 



(26) 
(27) 



where /z's are in units of i. In the case J ^ 0, the last 
term in Eq. (j2"3"|) near the ±(&d, fco)-nodes is expanded as 



^(4-/i)^|| 

and near the ±(fco, ~ fc,o)-nodes 

± Jq — vjSk±, 



(28) 



(29) 



where Jo = J(l — f )VW M 2 is the energy shift of the 
nodes, and = J(— ^i 2 + 2^i + 8). As a result of the 
shift, two Fermi pockets are induced, as shown in FigfTJ 
giving a finite density of states at zero energy. 



C. Continuum formulation and the linearized 
Hamiltonian 

In the low temperature specific heat measurement, 
only the low energy excitations contribute to the result. 
QP's near the ±(fc_D,fc_o) nodes may be expected to re- 
sult in the \f~H behavior of the Sommerfeld coefficient 
7(i?). On the other hand, the low energy QP's near the 
Fermi surfaces at ±(&d, — &d) may be expected to give 
rise to the finite zero field Sommerfeld coefficient 7(0) 
and perhaps even the quantum oscillations in high fields. 
To test this, we formulate the continuum version of the 
BdG Hamiltonian in the presence of loop current order, 
and linearize it near the four nodal points. 

In the absence of loop current order, the continuum 
Hamiltonian readsi 3 . 



u 



A* 



A 

- ft- 



(30) 



with Tie = l/2m(p — e/cA) 2 — /1, p = — iftV the mo- 
mentum operator, and V x A = Hz. In the following, 
we choose the a;-axis along the (kD,kn) nodes and the 
y-axis along the (—/cd, ko) nodes. Then the gauge in- 
variant d-wave pairing operator has the form 

A = 4-{px, {Py, A(r)}} + -VA(r)fep^), (31) 



Pf 



with pf the Fermi momentum and 4> the phase of the 
superconducting gap A(r). The curly bracket represents 
symmetrization, {a, b} = \/2(ab + ba). After the singular 
gauge transformation (fT2]) , the Hamiltonian become o 12 ' 13 



b 



where D = A /2p 2 F [p x 



D 



-2k(p- mv f ) 2 + / 1 



(32) 



a y ] + [x «->• y) and 
l/m(hV(f)^ — e/cA) for [i = A,B. The Berry vec- 
tor potential! 2 , a = to/2(v^ - vf ) = h/2(\7cj) A - V<p B )- 
The linearized approximation in the vicinity of one of the 
±(fc£>, — kr>) nodes results in 



H 



N 



"Ho + H j 



(33) 
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where 

Ha = ( VFf>y V ^ ) (34) 

\ V APx -VFPy J 

is the free Dirac Hamiltonian and 

n'=( mVF < VA \). (35) 

V "a«i mv F vf y J 

In the above vf is the Fermi velocity and va = ^-o/Pf is 
the slope of the gap at the node. Hn can be written as 

T-Ln = v F {p y + a y )a 3 + va{Px + a x )ai + mv F v sy , (36) 

where v s = (vf + vf )/2 = l/m(h/2\7(j) - e/cA) is the 
superfluid velocity. From Hn it is readily seen that a 
couples to the Dirac fermions as a vector potential while 
v s results in a Doppler shift. The magnetic field pro- 
duced by a consists of a set of ±7r-fmx delta function 
spikes at the vortex cores and vanishes on average. It 
does not lead to Landau level quantization^ 2 -. 

It is expected that physical quantities should be in- 
dependent of the choice of A and B sublattices, since 
there should be no physical distinction between A and 
B vortices. However, as discussed in Ref.—, two distinct 
choices of A-B sublattices as illustrated in Fig. 2 of Ref. 13 
result in qualitatively similar but still somewhat differ- 
ent band structures and densities of states, particularly 
at higher energies. Despite significant effort— this prob- 
lem remains a bit of a challenge: while the large gauge 
invariance is easily restored by judicious enforcement of 
boundary conditions at vortex locations, the interference 
among the nodes in a perfect vortex lattice obscures the 
ultimate choice for these boundary conditions^. At any 
rate, these mathematical subtleties are inherent only to 
the linearized BdG Hamiltonian and they do not arise at 
all in the tight-binding lattice formulation. 

The linearized Hamiltonian associated with loop cur- 
rent order can be derived from Eq. (fT9|) . Near one of the 
Fermi pockets, it reads 

Hj = {-vjp y + J )t. (37) 

In a magnetic field, after the singular gauge transforma- 
tion, it becomes 

- vj(p y + a y )t - mvjv S ya 3 + J t. (38) 

Therefore, the full linearized Hamiltonian near one of the 
Fermi pockets is 

H = vf(P v + a y )a 3 + va(Px + a,x)cri + mvFV S y 

- Vj(j)y + CLy) - TnVjV S y(T3 + Jq , (39) 

from which it is seen that the superfluid velocity couples 
to the Bogoliubov QP's, in part, as a vector potential 
through loop current order, 

Vj tt G 

(v F p y - mvjv sy )a 3 = v F {p y -(^5 y A y ))a 3 (40) 

vf ^ c 



The effective vector potential a e ff has a zero x- 
component, while the y-component is (vj/vF)(c/e)mv sy . 
The associated effective field, b e ff = V x a e fj , vanishes on 
average in the magnetic unit cell. The linearized Hamil- 
tonian near the ±(fco, fco)-nodes, which are not shifted 
by the loop order, resembles Eq. (|39|) . but with Jq = 0. 



III. NUMERICAL RESULTS 

At low temperature, the non-oscillatory part of the 
specific heat C(T,H) is linear in the temperature T, for 
a Fermi liquid composed of Schrodinger particles with 
p 2 /2m dispersion. The Sommerfeld coefficient can be 
defined as j(H) = C(T,H)/T. Experimentally, it has 
been found that the non-oscillatory part of j(H) goes 
as VH in low field, which is consistent with the d-wave 
vortex state scenario. In the high field j(H) also develops 
an oscillatory component, which obeys LK formula^ 

C sc(T, H) = 

-AT±R DM ^)co S (2M£- c -l))r(x) 
p=i L L 

(41) 

where A is a constant, Rp — exp ((— 2 / n 2 pkBT£,)/ (fujj c )) 
is the Dingle factor, lu c — eH/(m*c) is the cyclotron fre- 
quency with m* the effective mass, x — 2ir 2 pk bT '/ '(fiw c ), 
f"(x) = x( (l + cosh 2 x)j sinh 3 x — 2 cosh a;/ sinh 2 x), Jo is 
the Bcssel function of the first kind (not to be confused 
with the energy shift by loop current order), and t w is 
the c-axis hopping energy. In the experiments only the 
first harmonic with p — 1 is identified. In the presence of 
loop current order, J ^ 0, and the Sommerfeld coefficient 
becomes j(H,J). We will compare our results with the 
LK formula and show that in the vortex state with loop 
current order the formula does not hold. 

In the appendix, we derive the formula for the os- 
cillatory part of the specific heat assuming that Dirac 
particles with velocities vf and va couple minimally to 
the vector potential corresponding to a uniform magnetic 
field. As we stressed before, the (i-wave Dirac particles 
do not have such coupling. Nevertheless, we find it useful 
to contrast our numerical finding to this analytical for- 
mula. In this case, the expression for C osc is similar to 
Eq. (j4T|) . but there is an important difference. The effec- 
tive mass in the above formula is replaced by Ep / [vfVa)- 
As such, the amplitude of the oscillations also depends 
on the Fermi energy in addition to the temperature and 
the magnetic field. 

We use realistic values of physical quantities of YBCO 
as parameters in our Hamiltonian. The Fermi velocity 
is taken to be vf = 2.15 x I0 5 m/s, the lattice con- 
stant a — 0.38nm, the doping 15%, and the Dirac cone 
anisotropy a = 14. We first get /i = 0.297i from the dop- 
ing, and then derive the nearest neighbor hopping energy 
t = 0.132ey using Eg . (f2"6"l) . Within our method, we are 
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not able to sweep the magnetic field continuously as men- 
tioned in the introduction. Instead, we sweep the Fermi 
pocket area by varying the loop current order or Zeeman 
energy in a fixed magnetic field. From this viewpoint, 
Onsager relation readsii 



A(^ +1 ) - A(^) = A 



(42) 



where is ^th energy level when a magnetic field is ap- 
plied. This means that the period of oscillations, which is 
the difference between the areas enclosed by the orbits of 
adjacent energy levels in fc-space, equals the area of the 
magnetic Brillouin zone A Q = = = (f^) 2 . 

We study the specific heat in four different magnetic 
fields, with magnetic length l& — 60a, 40a, 36a and 28a. 
For a = 0.38nm this corresponds to field strengths 7.7T, 
17. 4T, 21. 5T and 35. 6T, respectively. Being fully aware 
of the caveat that for Dirac particles the amplitude of the 
oscillations of the specific heat may also depend on the 
Fermi pocket area, and therefore strictly speaking does 
not follow the Onsager relation, we investigate whether 
such relation holds in the d-wave superconducting state 
with loop current order. 

We use Arnoldi algorithm to diagonalize the Hamilto- 
nian (fTTj)) . Only the low energy bands need to be taken 
into account since the high energy bands give negligible 
contribution to the low temperature specific heat. Us- 
ing t — 1580K, all bands below 100K are considered. 
This gives us enough accuracy to determine the specific 
heat up to ~10K. We find that a 40 x 40 mesh in the 
first magnetic Brillouin zone (corresponding to a system 
with 40 x 40 magnetic unit cells) gives convergent results, 
showing little difference from that with a 80 x 80 mesh 
at the temperatures under study. 

In what follows, we present results for the Sommer- 
feld coefficient with loop current order j(H, J) in fixed 
magnetic fields H while J is continuously swept, which 
are later compared with the results from sweeping the 
Zeeman energy. 



A. Oscillations as a function of loop current order 

1. Frequency of oscillations 




0.02 



0.04 



0.06 



0.08 



J[t] 



FIG. 3: The dependence of Fermi pocket area Af on loop 
current order J in the d-wave superconducting state. The 
dots calculated from numerics are fitted by a parabola. Af is 
in unit of a -2 , where a = 0.38nm is the linear size of the unit 
cell, and J is in unit of the nearest neighbor hopping energy 
t. The anisotropy a = 14. The interval between adjacent 
horizontal grid lines is Aq — (2it/£b) 2 , where the magnetic 
length Ib = 40a in this figure. There are about 8 intervals 
between J = 0.044 and 0.07t, marked by the vertical lines. 
(Inset) The electron-like and hole-like Fermi pockets in the 
1st Brillouin zone at J = 0.05t. 
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FIG. 4: The Sommerfeld coefficient y(H,J) as a function 
of loop current order J in a d-wave superconductor, in four 
different finite fields with magnetic length Ib = 28a (red), 36a 
(purple), 40a (blue) and 60a (green) and zero field (black). All 
curves are at 2K except that the Ib = 60a curve is at IK. 



In the presence of loop current order J, two Fermi 
pockets appear in the {ko, —ko) direction of the Bril- 
louin zonei^, as shown in FigfT] At small J, the area 
of each Fermi pocket Ap is quadratic in J since the 
dispersion has a Dirac cone structure. FigJ3] shows 
numerically calculated Ap vs. J for the tight bind- 
ing model. The interval between two adjacent horizon- 
tal lines is A$ = (2tt/£b) 2 , with the magnetic length 
Ib = 40a. We vary the Fermi pocket area by choosing 
J = <V0.04 2 + 0.5 x 10~ 4 n, with the integer n ranging 
from 1 to 70. Then J ranges approximately from Q.OAt 
to 0.07i, within which there are about 8 intervals of A . 

FigfJ] shows j(H, J) vs. J at low temperature for the 
four different finite fields and zero field. All the curves 



are at 2K except the one in the lowest field (green) which 
is at IK. The zero field 7(0, J) is calculated using the dis- 
persion Eq.p3|. For this temperature, the frequency of 
the oscillations basically obeys Onsager relation Eq. (l4"2l , 
i.e. the frequency (or the number of periods) of oscilla- 
tions is proportional to the inverse of the magnetic field. 
A comparison of the oscillations at two different fields 
as a function of the rescaled Fermi pocket area is shown 
in Fig JSJ where the zero field background has been sub- 
tracted. There are about 1.4 periods between two ad- 
jacent vertical lines for both fields, suggesting that the 
Onsager relation holds approximately, albeit with an ef- 
fective charge e* w 0.7e. 

The experimental results show that the background on 








FIG. 5: The Sommerfeld coefficient y(H,J) as a function of 
the rescaled Fermi pocket area Af/Aq, for is = 40a (blue) 
and 28a (red) at 2K after the zero field background is sub- 
tracted. 



FIG. 7: The Sommerfeld coefficient y(H,J) as a function of 
loop current order J at temperature from 5K to IK with step 
IK, in the field with magnetic length Ib = 40a. The phase 
shift basically agrees with that predicted by LK formula. 




FIG. 6: The Sommerfeld coefficient y(H, J) as a function 
of loop current order J, for four different finite fields with 
magnetic length Is = 28a (red), 36a (purple), 40a (blue) and 
60a (green) and zero field (black) at temperature 5K. In the 
inset, the four dots correspond to the intersection of the four 
finite field curves and the vertical dashed line in the main 
figure. For comparison, the solid line in the inset shows \/H 
dependence; the units of the vertical axis are the same as in 
the main figure. 

top of which the oscillations reside has a \/~H behavior— 
In FigJS] we show y(H, J) vs. J for the four finite fields 
and zero field at 5K. Although increasing with H, the 
high temperature y(H, J) deviates from yH behavior, 
as shown in the inset of FigJU As J increases, the finite 
field curves become closer to each other, but farther from 
the zero field curve. 



2. Temperature dependence of the oscillations 

The temperature dependence of the quantum oscilla- 
tion of j(H, J) is shown in Figf7]for i B = 40a =$> H s» 
17. 4T. At T > 5K, no oscillations appear. As the temper- 
ature is lowered, the oscillations arise and the amplitude 
grows with decreasing temperature. A phase shift is ob- 
served at IK. In LK formula Eq. (|4"Tj) . a phase shift for 
p = 1 occurs at f"(x) = 0, where x = 2ir 2 kBT/(huj c ). 



Here uj c = eH/(m*c). If we had charged Dirac particles, 
m* should be replaced by Ef/(vfva)- If the oscillations 
obey LK formula, then with different parameter config- 
urations, the phase shift should occur at x « 1.6. Plug- 
ging in the parameters T — IK and H = 17. 4T, using 
Ep = J where J is defined in Eq. (f29|) , and e* = 0.7e, we 
find that the phase shift should be located at J w 0.05i. 
This basically agrees with the numerical result in Fig|7] 

At small J, x ~ JT/H. Therefore, for a fixed mag- 
netic field H, the phase shift should occur at the value 
of J proportional to l/T. Similarly at different magnetic 
fields H, if T is fixed, the phase shift should occur at the 
value of J proportional to H . However, the phase shift 
can not be well identified in all the cases, thus we are not 
able to accurately test whether this relation holds. 



3. Dependence of the oscillations on the configuration of 
the vortex lattice 

We also study the oscillations in a vortex lattice with 
the two vortices placed at different positions in a mag- 
netic unit cell, while the size of the magnetic unit cell is 
kept the same. We choose Ib — 40a. The comparison is 
shown in FigJSJ The blue curve is the result with vortices 
distributed uniformly, with the separation Ib/2 = 20a in 
both the horizontal and the vertical directions (see Fig [2]), 
the same as the blue one in FigH] the red curve is the re- 
sult with the two vortices in the same magnetic unit cell 
placed much closer, with the separation 4a in both di- 
rections. The same frequency of oscillations is observed, 
excluding the possibility that the oscillations are from 
Bragg plane reflections due to specific vortex configura- 
tions. There is a difference between the phases of the two 
configurations, which suggests that in the resistive state 
with creeping vortices, the phase difference may smear 
out the oscillations. 
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FIG. 8: Oscillations of Sommerfeld coefficient, for different 
configurations of the vortex lattice. The separation between 
two vortices within a magnetic unit cell in both directions is 
20a (blue) and 4a (red), in the field with £ B = 40a. 
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At zero temperature, the oscillations of the Sommer- 
feld coefficient simply corresponds to the oscillations of 
the density of states (DOS) at zero energy. We now inves- 
tigate whether the DOS resembles that of Dirac fermions 
minimally coupled to an external magnetic field H. This 
can help us to illustrate the difference between the two 
systems. The DOS as a function of Fermi pocket area 
in the field H with i B = 40a is shown in Fig® The 
mesh shows the DOS of Dirac fermions with the same 
anisotropy and Fermi velocity minimally coupled to the 
same field H. Clearly, to a large extent, the d-wave su- 
perconductor with loop current order in the vortex state 
does not resemble Dirac fermions minimally coupled to 
a magnetic field. Therefore, the LK formula, which is 
derived for systems with Landau level quantization, does 
not hold in this case. Nevertheless, an oscillatory feature 
of the DOS is seen. It is this feature which is responsible 
for the oscillations of the specific heat at intermediate 
temperatures presented in Figs. (4-8). 

B. Oscillations as a function of Zeeman energy 

Up to now, we have neglected the Zeeman splitting 
due to the external magnetic field, since it only affects 
the results in an insignificant way and the conclusions 
do not change. We imagine changing the Zeeman term 
while holding the magnetic field and the pairing term 
fixed. This will also induce Fermi pockets in the d-wave 
superconducting state, even in the absence of the loop 
current order. In zero field, the dispersion with such a 
term is 

E k = ±y/t*+Al + E z (43) 

where the Zeeman term Ez shifts the energies, resulting 
in four Fermi pockets(see Fig[TU]). Since at small energies 
the pockets are ellipses, the area of each one Ap is easily 



FIG. 9: The DOS as a function of energy and the Fermi 
pocket area, A F , rescaled by Ao = *g-. Here £ B = 40a. The 
mesh shows the DOS of Landau quantized anisotropic Dirac 
fermions in the same field. The scale shown is such that the 
integrated DOS over the range of data is normalized to 1. 

calculated, and varies quadratically with Zeeman energy 
as shown in FigHOl We use the same parameters as in 
the loop current order case, and sweep Ez from 0.02£ to 
0.16£, for two magnetic fields with £g=28a and £s=40a. 
In FigQj] we show j(H) vs. Af/Aq. If Landau lev- 
els are formed and Onsager relation holds, the frequency 
of the two oscillations should be the same, regardless of 
the magnitude of the magnetic field. Nevertheless, the 
frequency is doubled when the field is doubled, which 
is consistent with QP's forming Bloch bands instead of 
Landau levels^. Comparing this result with the oscilla- 
tions induced by the loop current order, we conclude that 
the latter has a more intricate nature whose effects can 
not be accounted merely by the presence of the Fermi 
pockets. Rather, we believe, the special nature of the 
coupling between the loop current order and the nodal 
d x 2_ y 2 QP's is essential to account for the observed os- 
cillations obeying Onsager relation. 

IV. THE LINEARIZED PROBLEM 

In the continuum formulation we diagonalize the lin- 
earized Hamiltonian Eq. (|39|) in the plane wave basi s 12 ' 16 , 
and repeat the calculations above for two values of 
the magnetic field. The Fermi pocket area Ap — 
ttE f /(v f va), where E F = J given below Eq. (j2"fj|) . We 
find that the two pockets give the same contribution to 
oscillations of the specific heat, while the nodes do not 
contribute to the oscillations. Fig[T3] shows the oscilla- 
tions of j(H,J) in fields with Ib = 28a and Ib = 40a 
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FIG. 10: The dependence of Fermi pockets area Af on Zee- 
man energy in the d-wave superconducting state. (Inset) The 
four Fermi pockets at Ez = 0.02t (black) and Ez = 0.16t 
(red). The parameters are the same as in FigO 



FIG. 12: The Sommerfeld coefficient *y(H, J) as a function 
of J in fields with £ B = 28a(red) and £b = 40a (blue) at 5K 
resulting from one Fermi pocket in the linearized formulation. 
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FIG. 11: The Sommerfeld coefficient y(H) as a function of the 
rescaled area the Fermi pockets induced by Zeeman energy, in 
the magnetic fields with £ B — 40a (blue) and £b — 28a (red). 



from one of the pockets. With the same range of J as 
in Sec. Ill, the same frequency is found for both curves, 
which confirms our findings above. The effective charge 
e* e here. 

The linearized formulation enables us to sweep the ex- 
ternal magnetic field continuously at a fixed J, and to 
verify Onsager relation, which reads^i 



A(1) = ^L 

K H' hcA F 



(44) 



where A(l/H) is the period of oscillations if Ap is fixed 
and H is varied. This is equivalent to A{l 2 B ) = (27r) 2 /Ap. 
We take J = 0.05i and sweep £ B , the result of which 
is shown in FigUJJl The period is A£ 2 B w 200. Using 
Eq. ([26)) . (|44|) and the expression of Jo, we determine that 
the effective charge e* « e, which agrees with what is 
derived above for the tight-binding formulation. 



V. DISCUSSIONS AND CONCLUSIONS 

We have shown that the Fermi pockets induced by 
loop current order can give rise to quantum oscillations 



FIG. 13: 7(77, J) as a function of £% at J = 0.05t and T =2K, 
resulting from one Fermi pocket in the linearized formulation. 



of the specific heat within a limited temperature range 
which seem to obey Onsager relation, with an effective 
charge comparable to the electric charge. In the the- 
ory of metals, Onsager relation is established from argu- 
ments that in a uniform magnetic field, semiclassically, 
electrons move in constant energy surfaces with quan- 
tized energies^. In a d-wave superconductor with loop 
current order, however, such a argument does not work. 
Suppose that Bogoliubov QP's circulate the Fermi sur- 
face which is an ellipse. Although the Fermi surface is 
still a constant energy surface, the charge of a Bogoliubov 
QP varies with its position as it moves around the Fermi 
surface, and the average charge over the Fermi surface 
is zero^, which makes the argument for metals invalid 
here. 

We trace the origin of the oscillations observed in our 
numerical calculation to the highly inhomogeneous fic- 
titious magnetic field experienced by the d-wave Dirac 
particles in the presence of the loop order. Such coupling 
is absent if the Bogoliubov QP Fermi pockets are due 
to Zeeman effect only. Indeed in the temperature range 
where the oscillations appear, the Onsager relation holds 
approximately in the case with the loop order but does 
not hold if only the Zeeman shift is included. 

One disadvantage of this picture is that the contribu- 
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tion from the loop current order induced Fermi pockets is 
at odds with the experimentally confirmed \f~H behavior 
of the background on top of which the oscillations occur. 
Another disadvantage is that the full Lifshitz-Kosevich 
relation ceases to hold since ultimately there is no Lan- 
dau quantization as shown in Fig.(|9|). 

Recently, a preprint Ref.— with a similar concern ap- 
peared online. Comparing with their work, we agree 
that there are no quantized Landau levels near the Fermi 
pockets of Bogoliubov QP's. Nevertheless, we do find 
quantum oscillations of the specific heat in a small tem- 
perature range, while, by analyzing the density of states 
at the Fermi level, they conclude that no quantum os- 
cillations should appear unless spin stripe order coexists 
with loop current order. 



Acknowledgments 

We wish to thank S. A. Kivelson and C. M. Varma for 
useful discussions. This work was supported by a NSF 
CAREER award under Grant No. DMR-0955561, the 
NSF Cooperative Agreement No. DMR-0654118 and by 
the State of Florida . 



then 



[a,a T ] = 1 



and the Hamiltonian can be written as 
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The eigenenergies of this Hamiltonian are 
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where C = V% J^ig - We consider only the positive 

eigenenergies e„ = (,\fn which are connected to the neg- 
ative eigenenergies by particle-hole symmetry. 
The grand potential is then 



fl = -TD^ln(l +, 
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Appendix A: Quantum oscillations of the specific 
heat of Dirac fermions 



where D = gHL 2 /cf>o is the degeneracy of one Landau 
level, g the number of species of fermions. Since 



In this appendix we consider a hypothetical problem, 
the two-dimensional anisotropic Dirac fermions moving 
in a perpendicular magnetic field Hz. We use it to set up 
a contrast with the results found for the physical system 
considered in the main text. The Hamiltonian reads 



then 



n=o n=0 ^ * 
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and olf) is the anisotropy. Since 
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(Al) 



(A2) 
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where C-b — y jjj (defined differently from that in the 
main text), if we define the annihilation and creation 
operator as 



,t - 
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y/2n olf 
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Integrating by parts, we have 
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The first two terms are non-oscillatory. Integrating by Using the formula 
parts again for the third (oscillatory) term, we get 
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and differentiating the grand potential twice with respect 
to T, we arrive at the oscillatory part of the Sommerfeld 
(A19) coefficient 
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(A21) 



The first exponential accounts for the quantum oscilla- 
tions of the Sommerfeld coefficient with the Fermi pocket 
area and the magnetic field, and the integral determines 
the amplitude as well as the phase shift of the oscillations. 

At the low-temperature limit, fi/T 3> 1, we can extend 
the lower limit of the first integral to — oo since the dif- 
ference is exponentially small. Also, in the integral, the 
first exponential oscillates much faster than the second 
one since x is bounded by 1 / cosh 2 | to <~ 1 . Thus we 
drop the second exponential as well as the last factor, 
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The integral has the analytical form 
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FIG. 14: The integral appearing in Eq. (|A22|) 
(blue) and Eq. (|A"2il (red), with T/£ = 0.1. 
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The Fermi pocket area is 
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So the combinations appearing in the Sommerfeld 
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The coefficient is exactly the same as in the Schrodinger 
case, except that m* has a different definition and de- 
pends on \i now. In the Schroding er case, C osc /T is 
strictly periodic in the Fermi pocket area since the am- 
plitude has no dependence on /x, but for Dirac fermions 
it is different. Because z depends on ^i, the amplitude, 
as well as the phase shift, depend not only on the tem- 
perature and the magnetic field, but also on the Fermi 
pocket area. 

From the first harmonic p = 1, we can deduce the 
Onsager relation: the period of the cosine as a function 
of Ap and 1/H, respectively, is 



A(A F ) 



Att 2 H 



4n 2 



(A34) 
(A35) 



To compare the approximate expression of the amplitude 
obtained in the limit fi/T 1 with the expression valid 
to any fi/T, in FigHjwe plot the integral in Eq. (|A22p 
and that in Eq. ([A21[ ) as a function of n/T, while fixing 
T/( — 0.1. The two curves basically coincide at /i/T > 
10. The position where the phase shift occurs, i.e. where 
the integral changes sign, is different for the two curves. 
At a smaller T/£, the phase shift occurs at a larger [ifT 
for both curves, and the positions where the phase shift 
occurs become closer. 
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